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METHOD FOR ANALYTICALLY OBTAINING CLOSED FORM EXPRESSIONS 
FOR SUBSURFACE TEMPERATURE DEPTH DISTRIBUTION ALONG WITH ITS 
ERROR BOUNDS 

Field of the invention 

The present invention relates to analytical solution to error bounds on the subsurface 
temperature depth distribution. This invention has a wide range of application in quantifying 
the subsurface thermal state of the crust and has a wide range of applications. The method of 
the invention has helped in obtaining closed form expressions for the subsurface temperature 
depth distribution along with its error bounds. The exact formulae are useful in better 
evaluating the thermal state and have a wide range of applications in oil and natural gas 
prospecting, tectonic studies and mineral prospecting. 
Background of the invention 

This invention relates for obtaining and computing the subsurface temperature depth 
distribution along with its error bounds. The solution has been determined for the stochastic 
heat conduction equation by considering different sets of boundary conditions and radiogenic 
heat sources and incorporating randomness in the thermal conductivity. In understanding the 
Earth thermal structure there are several questions which need clear answers. Many of the 
controlling parameters that define the Earth's processes are not known with certainty. In such 
situations these controlling parameters can be defined in a stochastic framework and an 
average picture of the system behavior together with its error bounds can be quantified. 

The thermal structure of the Earth's crust is influenced by its geothermal parameters 
such as thermal conductivity, radiogenic heat sources and initial and boundary conditions. 
Basically two approaches of modeling are commonly used for the estimation of the 
subsurface temperature field. These are: (1) deterministic approach and (2) the stochastic 
approach. In the deterministic approach the subsurface temperature field is obtained assuming 
that the controlling thermal parameters are known with certainty. However, due to 
inhomogeneous nature of the Earth's interior some amount of uncertainty in the estimation of 
the geothermal parameters are bound to exist. Uncertainties in these parameters may arise 
from the inaccuracy of measurements or lack of information about the parameters 
themselves. Such uncertainties in parameters are incorporated in the stochastic approach and 
an average picture of the thermal field along with its associated error bounds are determined. 
To assess the properties of the system at a glance we need to obtain the mean value that gives 
the average picture and the variance or the standard deviation that is the variability indicator 
which gives the errors associated with the system behavior due to errors in the system input. 
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Subsurface temperatures are also seen to be very sensitive to perturbations in the input 
thermal parameters and hence several studies have been carried out in quantifying the 
perturbations in the temperatures and heat flow using stochastic analytical and random 
simulation techniques. Quantification of uncertainty in the heat flow using a least squares 
inversion technique incorporating uncertainties in the temperature and thermal conductivities 
has been done, Tectonophysics, Vol 121, 1985 by . Vausser et al . The effect of variation in 
heat source on the surface heat flow has also been studied, Journal Geophysical Research, V 
91, 1986, by Vasseur and Singh, Geophysical Research Letters, V14, 1987, by Nielsen. In 
most of the studies the stochastic heat equation has been solved using the small perturbation 
method. Using the small perturbation method the heat conduction equation has been solved 
by incorporating uncertainties in the heat sources and the mean temperature field along with 
its error bounds have been obtained, Geophysical Journal International, 135, 1998, by 
Srivastava and Singh. The random simulation method has also been used to model the 
thermal structure incorporating uncertainties in the controlling thermal parameters, 
Tectonophysics, V156, 1988 by Royer and Danis, Marine and Petroleum Geology, V 14, 
1997, by Gallagher et al, Tectonophysics, V 306, 1999a,b, by Jokinen and Kukkonen. This 
numerical modeling is very useful in studying the nonlinear problems but sometimes simple 
1-D analytical solution to the mean behavior and its associated error bounds is very useful in 
quantifying the uncertainty. The stochastic differential equations in other fields are now being 
solved by yet another approach called the decomposition method, Journal of Hydrology, V 
169, 1995, by Serrano. In a recent study using this new approach the stochastic heat equation 
has been solved incorporating uncertainties in the thermal conductivity where the solution to 
the temperature field is obtained using a series expansion method, Geophysical Journal 
International, V 138, 1999, by Srivastava and Singh. The thermal conductivity is considered 
to be a random parameter with a known Gaussian colored noise correlation structure. 

In this invention the stochastic solution to the mean and variance in the temperature 
field for a different set of boundary conditions and different radiogenic heat source function 
has been obtained following the procedure of Geophysical Journal International, V 138, 
1999, by Srivastava and Singh. The expressions for mean and variance in temperature depth 
distribution for different heat sources and boundary conditions have been obtained and used 
to compute and plot the subsurface thermal field along with its error bounds. 
Objects of the invention 

The main object of the invention is to provide analytical solution to error bounds on 
the subsurface temperature depth distribution which obviates the drawbacks detailed above 
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Another object of the invention is to provide an efficient method for obtaining closed 
for solution to error bounds on temperature depth distribution for different set of boundary 
conditions. 

Yet another object of the invention is to provide for quantification of subsurface 
temperature depth distribution and its error bounds for known Gaussian thermal conductivity 
structure. 

Brief Descriptions of the Drawings 

The attached figures ( Fig 2-7) are the plots of mean temperature ± 1 standard 
deviation for different values of the controlling thermal parameters. 
Fig. 1 is a flow sheet depicting the method of the invention. 

Fig 2 Plot of mean temperature ± 1 standard deviation for condition 1 when no heat 
source is considered and the boundary conditions employed are constant surface temperature 
and constant surface heat flow. 

Fig 3 Plot of mean temperature ± 1 standard deviation for condition 2 when no heat 
source is considered and the boundary conditions employed are constant surface temperature 
and constant basal heat flow. 

Fig 4 Plot of mean temperature ± 1 standard deviation for condition 3 when a 
constant heat source is considered and the boundary conditions employed are constant 
surface temperature and constant surface heat flow. 

Fig 5 Plot of mean temperature ± 1 standard deviation for condition 4 when a 
constant heat source is considered and the boundary conditions employed are constant 
surface temperature and constant basal heat flow. 

Fig 6 Plot of mean temperature ± 1 standard deviation for condition 5 when an 
exponential heat source is considered and the boundary conditions employed are constant 
surface temperature and constant surface heat flow. 

Fig 7 Plot of mean temperature ± 1 standard deviation for condition 6 when an 
exponential heat source is considered and the boundary conditions employed are constant 
surface temperature and constant basal heat flow. 
Summary of the invention 

Accordingly the invention provides an analytical solution to error bounds on the 
subsurface temperature depth distribution, which comprises a method of solving the heat 
conduction equation incorporating Gaussian uncertainties in the thermal conductivity. 
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In an embodiment of the present invention the stochastic heat conduction equation has 
been solved using a series expansion method to obtain the closed form solution to the mean 
and variance in the temperature depth distribution. Simple deterministic solution to the 
problem is not sufficient and quantifying the errors in the system output due to errors in the 
input parameters is very essential. These errors bounds are very important for a better 
evaluation of the subsurface thermal structure. 

In another embodiment of the present invention the expression for mean temperature 
and the variance in temperature are obtained for six different set of prescribed boundary 
conditions. 

Detailed description of the invention 

The present invention deals with the solution to the governing stochastic heat 
conduction equation to obtain the mean and variance in the temperature fields as shown in the 
flow chart. The heat conduction equation with random thermal conductivity is expressed as 

±{K + K'(z))^f- } = -Mz) (1) 
dz dz 

where 

T is the temperature ( °C) , 

A(z) is the radiogenic heat source ( {iW / m 3 ) , 

K(z) = K + K (z) is the thermal conductivity ( W I m °C) 

which is expressed as a sum of a deterministic component and a random component 

K (z) is the random component with mean zero and a Gaussian colored noise correlation 
structure represented by 

E(K(z)) =0 (2) 

E(K( Zl )K(z 2 )) = a 2 K e- p|2l - Z2 > (3) 
where 

a 2 K is the variance is thermal conductivity ( W I m °C) 2 

p is the correlation decay parameter m ~ l (or 1/p is the correlation length scale) 

zi and Z2 are the depths (m) 

Following the procedure of given in Geophysical J International, V 138, 1999 by 
Srivastava and Singh, the solution to mean temperature and its standard deviation has been 
obtained for three conditions of heat sources (1) Zero (A(z) = 0) (2) Constant (A(z) = A) and 
(3) Exponentially decreasing with depth ( A(z) = Age'^ 0 ) and 
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Associated boundary conditions are defined by 
Type(i) Boundary Condition : 
Constant surface temperature 

T=T 0 at z=0 (4) 

Surface heat flow Q 5 ( mW / m 2 ) 

K^ = Q S at z=0 (5) 
dz 

Type (ii) Boundary Condition : 
Constant surface temperature 

T = T 0 at z = 0 (6) 
Basal heat flow Q B (mW I m 2 ) 

K^f = Q B at z = L (7) 

dz 

The standard deviation, which is a measure of error in the system output, is obtained 
by taking the square root of the variance. The solutions to different conditions obtained have 
been given below. 

Condition 1: When no heat source is considered and the B.C used are the surface 
temperature and surface heat flow 

The governing heat conduction equation without heat source term is 

JL{ ( k + K(z)) *L } =0 (8) 
dz dz 

with constant surface temperature and constant surface heat flow as boundary conditions 

(Type (i)) 

Solution 

Mean Temperature 

f = E( T(z) )= T 0 +Q,z (9) 

K 

Variance in temperature 

a 2 = cl xTerml 

(10) 

where 

cl= 2C\Q) /K 2 
and 
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Terml = p^- - — + -^{l - e'<*) - 

where C K = a K / K is the coefficient of variability in the thermal conductivity. 
Condition 2: When no heat source is considered and the B.C used are surface 
temperature and heat flow at the base of the model ( Q s = Q B ) 

The governing heat conduction equation without heat source term is 

±{(K + K'(z))^- } =0 (11) 
dz dz 

with constant surface temperature and constant basal heat flow as boundary conditions 

(Type (ii)) 

Solution 

Mean Temperature 

f = E(T(z))= T 0 ^^-z (12) 

Variance in temperature 

<j 2 T =clxTerml (13) 

where 

cl= 2C 2 K Q B 2 IK 2 
and 

Term 1 is same as given in condition 1. 

Condition 3: When constant heat source is considered and the B.C used are the surface 
temperature and surface heat flow 

The governing heat conduction equation with constant heat source term is 

JL{ ( k + K (z) )*L } = _ A (14) 
dz dz 

with constant surface temperature and constant surface heat flow as boundary conditions 

(Type (i) ) 

Solution 

Mean Temperature 

f = E( T(z))= T 0 +%z ~ 4|r (15) 

Variance in temperature 

a 2 . =clxTemi +c2xTen?2 + c3xTerrrt + c4xTen?4 (16) 
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where 

cl = C\(A-pQ s f /K 2 
c2 = C 2 K Ap( A-pQ s )/K 2 
c3 =c2 

c4 = C\ p 2 A 2 / K 2 
and 

Terml = -L ( p^- + L1(- & - 1 + e ") - ( -^^-( pz + e* - 1 ) ) 
P 3 p 2 p 2 

Term! = ( -^±( (( pz + e~" - 1 ) / p 2 ) + p^ +4 " ^ " *' > ' 



p J 2 2 3 12 p 

(p^_2p)zi_ ^ + <Pz ~2)z 2 _ (JX + 2) -„ + 

p 3 6 12 p 2p 3 p 5 



Tenn3 = "Vf ^- + « /* + ' I*"* -1)^-^ ) - pz 2 e~>* ) 



P P P 

Tenn4 = (£±^L ( ^L_^. + uV -)/p + ((- (yQz + i) e - + i)(i^)--4.) 

/> 2 12 6 yp 2 20/? 

P /> P P 

12 6 /? 3 20/? /? 2 

Condition 4: When constant heat source is considered and the B.C used are the surface 

temperature and heat flow at the base of the model 

The governing heat conduction equation with constant heat source term is 

^-{(K + K'(z) } =-A (17) 

dz dz 

with constant surface temperature and constant basal heat flow as boundary conditions 

(Type (ii)) 

Solution 

Mean Temperature 

f = E(T(z))= T 0 + (Qb+ ** L) z -44 ( 18 > 

K. 2/C 

Variance in temperature 

=clxTerni + c2xTerr?S + c3xTern3 + c4xTern4 (19) 
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where 

cl= C 2 K iA-p(Q B +A*L)) 2 IK 2 
c2= C 2 K Ap(A- p(Q B +A*L))I K 2 
c3 =c2 

c4= C 2 K p 2 A 2 / K 2 

Terml , Term2 , Term3 and Term4 are same as given in condition3. 

Condition 5: When an exponential heat source function is considered and the B.C used 
are the surface temperature and surface heat flow 

The governing heat conduction equation with an exponential heat source term 

±{ ( K + K(z))lf- } = -vr 2/D (20) 

dz dz 

with constant surface temperature and constant surface heat flow as boundary conditions 

(Type (i)) 

Solution 

Mean Temperature 



f = E(T(z))= T 0+ 2±z + ^(l--^-^ /D j 



(21) 



Variance in the temperature 

crj. =clxTerni + c2xTern2 + c3xTerml + c4xTern4 (22) 
where the constants are 
cl= C\ AlO—pD) 2 I K 2 
c2 = C\ \p (pD - 1)( Q s - \D) I K 2 
c3 =c2 

c4 = C\p 2 (Q s -A 0 D) 2 IK 2 

The closed form solution for the integrals the in the above equation are 

Terml = - — {( pD — 1 )( 2z 2 — 2zD — D 2 e~ 2z/D + D 2 ) 

4(p-l/D) 2 C 



4[z(p-l/D) + l) r _ A , 7 m _ m -^ im 
(P + UD) 



2 [-z(p + l/D)-e- z ' p+uu > +1] 



+ [2zD+ D 2 e' 2zlD - D 2 ] } 

ApD + \)(2z 2 -2zD-D 2 e- 2z,D + D 2 ) 

4(yO + l/D) 2 1 
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+ -[-z{p-HD) e -W» +e -™° - e ~WO) j 

(p-l/D) 2 



- [ 2zD + D 2 e' 21 " 3 -D 2 ] } 
Term2 = \ { 2p ( z 2 D - 2zD 2 - 2D 3 e' z/D + 2D 3 ) 



p 2 



(l + PZ \ [z(p + l/D) + e-^ +UD) -1 ] 
(p + l/D) 2 



(p-voy 



[-z(p-l/D) + e zip ' UD) -1 ] } 



Termh = T {(p - l/D)( Z 2 D - 2zD 2 - 2D 3 e' z/D + 2D 3 ) 

(p-l/D) 2 1 

- Z ^- l/ 2 D) + l lpz + e-*-U 
P 

+ [zD+ D 2 e' z/D -D 2 ] } 

+ A(p + ll D)(z 2 D-2zD 2 - 2D 3 e~ tlD + 2D 3 ) 

(p + l/D) 2 lH 

-z(.p+l/D) 

[-pz + e<* -1 



P 2 



- [ zD+ D 2 e- z/D -D 2 ] } 



TermA = \ 2 3_(^+0 
P [3 P 



e 



P 2 



[-pz + e<* -I] } 



Condition 6: When an exponential heat source function is considered and the B.C used 
are the surface temperature and heat flow at the base of the model 

The governing heat conduction equation with an exponential heat source term 

±{(K + K(z))^- } = -\>e- z,D (23) 
dz dz 

with constant surface temperature and constant basal heat flow as boundary conditions 

(Type (ii)) 

Solution 

Mean Temperature 



T = E(T(z)) = T 0 +Si-z + ^£-(\-^ e - LID -e- z,D ] 
0 K K \ D J 



(24) 
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Variance in the temperature 

o$ =clxTemi +c2xTern2 + c3xTern3 + c4xTerr>4 (25) 
where the constants are 
cl = C\ A%(l-pD) 2 I K 2 

c2= C 2 K AqP (pD — 1)( Q g — A^De~ LfD ) I K 2 
c3 =c2 

c4= C 2 K p 2 (Q B -\De- L!D ) 2 IK 2 

The terms Terml, Term2, Term3, Term4 are same as given in condition 5. 
The present invention has its novelty over previous work in the following counts 

1. The method uses randomness in the thermal conductivity structure to quantify the errors 
in the subsurface temperature depth distribution 

2. The method used has led to exact closed form solution to the mean and its variance on the 
subsurface temperature field for different prescribed boundary conditions. 

3. The present solutions will be used extensively for quantifying the subsurface 
temperatures for any given region.. The exact formulae for mean and variance in the 
subsurface temperature depth distribution have not been given so far and have a wide 
application in geothermal studies. 

The following examples are given by way of illustrations and therefore should not be 
constructed to limit the scope of the present invention. 
Example 1 

Numerical values of the controlling input thermal parameters for a Realistic Earth 

model for condition 1 

Boundary conditions : 

Surface Temperature ( T 0 ) 0 (° C ) 

Surface heat flow ( Q s ) 80 ( mW / m 2 ) 

Model Depth (L). 3.5 (km) 

Random thermal conductivity : 

Mean thermal conductivity K 3.0 ( mW/m 2 ) 

Coefficient of variability C k 0.4 
Correlation length scale 1 / p 1 km 

Using these controlling thermal parameters the mean temperature and its error bounds 
have been computed using equations(9) and (10) and the results have been plotted in figure 1. 
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From the figures we see that the error bounds on the temperature increases with depth, 
increases with an increase in the coefficient of variability in thermal conductivity and in the 
correlation length scale 
Example 2 

Numerical values of the controlling input thermal parameters for a Realistic Earth 
model for condition 2 
Boundary conditions : 

Surface Temperature ( T 0 ) 0 (°C ) 

Basal heat flow ( Q B ) 30 (mW/m 2 ) 

Model Depth fL) . 5.5 (km) 

Random thermal conductivity : 

Mean thermal conductivity K 2.0 ( mW/m 2 ) 

Coefficient of variability Ck 0.3 

Correlation length scale 1 / p 1 km 

Using these controlling thermal parameters the mean temperature and its error bounds 
have been computed using equations(12) and (13) and the results have been plotted in figure 
2. From the figures we see that the error bounds on the temperature increases with depth, 
increases with an increase in the coefficient of variability in thermal conductivity and in the 
correlation length scale 
Example 3 

Numerical values of the controlling input thermal parameters for a Realistic Earth 

model for condition 3 

Boundary conditions : 

Surface Temperature ( T 0 ) 30 ( ° C ) 

Surface heat flow ( Q s ) 40 ( mW/m 2 ) 

Heat Source 

Radiogenic heat production (A) 2.5 ( juW/m 3 ) 

Model Depth (L). 10 (km) 

Random thermal conductivity : 

Mean thermal conductivity K 3.0 ( mW/m 2 ) 

Coefficient of variability Ck 0.2 
Correlation length scale 1 / p 3 km 
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Using these controlling thermal parameters the mean temperature and its error bounds 
have been computed using equations(15) and (16) and the results have been plotted in figure 
3. From the figures we see that the error bounds on the temperature increases with depth, 
increases with an increase in the coefficient of variability in thermal conductivity and in the 
correlation length scale. 
Example 4 

Numerical values of the controlling input thermal parameters for a Realistic Earth 

model for condition 4 

Boundary conditions : 

Surface Temperature ( T 0 ) 30 (°C) 

Surface heat flow ( Q B ) 20 ( mW / m 2 ) 

Heat Source 

Radiogenic heat production (A) 2.5 ( juW/m 3 ) 

Model Depth (L). 10 (km) 

Random thermal conductivity : 

Mean thermal conductivity K 3. ( mW/m 2 ) 

Coefficient of variability Cr 0.5 

Correlation length scale | | 4 km 

Using these controlling thermal parameters the mean temperature and its error bounds have 
been computed using equations(18) and (19) and the results have been plotted in figure 4. 
From the figures we see that the error bounds on the temperature increases with depth, 
increases with an increase in the coefficient of variability in thermal conductivity and in the 
correlation length scale. 
Example 5 

Numerical values of the controlling input thermal parameters for a Realistic Earth 
model for condition 5 

Boundary conditions : 

Surface Temperature (T 0 ) 30 (°C) 

Surface heat flow ( Q s ) 43 ( mW/m 2 ) 

Heat Source 

Radiogenic heat production ( A) 2.6 ( jjW /m 3 ) 

Characteristic Depth (D) 12 (km) 

Model Depth (L). 35 (km) 
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Random thermal conductivity : 

Mean thermal conductivity K 3.0 ( mW/m 2 ) 

Coefficient of variability Ck 0.3 

Correlation length scale 1 / p 10 km 

Using these controlling thermal parameters the mean temperature and its error bounds 
have been computed using equations(21) and (22) and the results have been plotted in figure 
5. From the figures we see that the error bounds on the temperature increases with depth, 
increases with an increase in the coefficient of variability in thermal conductivity and in the 
correlation length scale. 
Example 6 

Numerical values of the controlling input thermal parameters for a Realistic Earth 
model for condition 6 

Boundary conditions : 

Surface Temperature ( To ) 0 ( 0 C ) 

Surface heat flow ( Q B ) 20 ( mW/m 2 ) 

Heat Source 

Radiogenic heat production (A) 2.2 ( juW/m 3 ) 

Characteristic Depth (D) 10 (km) 

Model Depth (L). 35 (km) 

Random thermal conductivity : 

Mean thermal conductivity K 2.6 ( mW/m 2 ) 

Coefficient of variability C k 0.2 

Correlation length scale 1/ p 11 km 

Using these controlling thermal parameters the mean temperature and its error bounds have 
been computed using equations(24) and (25) and the results have been plotted in figure 6. 
From the figures we see that the error bounds on the temperature increases with depth, 
increases with an increase in the coefficient of variability in thermal conductivity and in the 
correlation length scale. 

THE MAIN ADVANTAGES OF THE INVENTION ARE: 

1. The advantage of this invention is that exact formulaes have been given to quantify 
the error bounds on the subsurface temperatures due to errors in the thermal 
conductivity for a conductive earth model. 
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2. The errors in the temperatures will help in a better evaluation of the crustal thermal 
structure. 

3. This study can be used in quantifying the conductive thermal structure along with its 
error bounds for any given region and thereby help in understanding the geodynamics 
of the region. 

4. The exact formulae for the mean temperature and its standard deviation can be used in 
a better evaluation of the thermal state of oil bearing regions. These exact solutions 
can be used in tectonic studies and in studies related to crystallization of minerals. 
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